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We describe collective-move Monte Carlo algorithms designed to approximate 
the overdamped dynamics of self-assembling nanoscale components equipped with 
strong, short-ranged and anisotropic interactions. Conventional Monte Carlo sim- 
ulations comprise sequential moves of single particles, proposed and accepted so as 
to satisfy detailed balance. Under certain circumstances such simulations provide 
an approximation of overdamped dynamics, but the accuracy of this approximation 
can be poor if e.g. particle-particle interactions vary strongly with distance or an- 
gle. The twin requirements of simulation efficiency (trial moves of appreciable scale 
are needed to ensure reasonable sampling) and dynamical fidelity (true in the limit 
of vanishingly small trial moves) then become irreconcilable. As a result, single- 
particle moves can underrepresent important collective modes of relaxation, such as 
self-diffusion of particle clusters. However, one way of using Monte Carlo simulation 
to mimic real collective modes of motion, retaining the ability to make trial moves 
of reasonable scale, is to make explicit moves of collections of particles. We will out- 
line ways of doing so by iteratively linking particles to their environment. Linking 
criteria can be static, conditioned upon properties of the current state of a system, 
or dynamic, conditioned upon energy changes resulting from trial virtual moves of 
particles. We argue that the latter protocol is better-suited to approximating real 
dynamics. 

I. INTRODUCTION 

Pairwise-interacting components in implicit solvent are often used as models of self- 
assembling systems, such as crystal-forming [TH3] and capsid-forming- [H [5] proteins, and 
patchy nanoparticles j6j [7] . In this paper we shall summarize ways of using Monte Carlo 
(MC) simulation to evolve such components in order to approximate the (interacting) Brow- 
nian motion that their real counterparts execute. Brownian motion is usually approximated 
in simulations by integration of overdamped equations of motion, called the Brownian Dy- 
namics (BD) method [HI E], with Monte Carlo algorithms more often used as a means of 
sampling thermal distributions. However, recent work shows that in certain circumstances 
the MC method can also evolve components according to an approximately correct dynam- 
ics. Dynamic MC simulations even offer some advantages over their BD counterparts: it 
is easier and computationally cheaper to evaluate potentials (MC) than forces (BD); MC 
can cope with pathological potentials (e.g. hard particles, abrupt changes in potential); and 
one does not face problems of numerical instability with MC as one does with BD (even for 
smooth potentials), and so can make larger basic moves. 



2 



In what follows we outline the reasons why single-particle Metropolis MC can effect a 
dynamics close to Brownian motion, and we summarize the work of others that makes use 
of this correspondence. We then describe extensions of the MC scheme that incorporate 
explicit moves of collections of particles, and argue that such schemes can be used to pre- 
serve the approximate realism of the MC method when collective modes of motion become 
important to a system's evolution. We take the view that although MC methods are not 
dynamically realistic in all details, they offer for some applications a convenient alternative 
to conventional integration of equations of motion. 



II. THE DYNAMICAL CHARACTER OF MONTE CARLO MOTION 

It is well-known that sequential moves of single components, proposed in an unbiased 
fashion and accepted according to the Metropolis criterion, permits sampling of the Boltz- 
mann distribution given sufficiently long simulation times [8j. However, is also true that if 
we restrict such moves to local translations and rotations then the dynamics executed by a 
single particle in an external forcefield is equivalent, in the limit of small trial moves, to a 
Langevin dynamics p^Ol [TTj . i.e. to Brownian motion in that potential. To see this, consider 
a particle in one dimension in a position-dependent potential U(x). The master equation 
corresponding to a Metropolis MC algorithm in which particle displacements x = x f — x are 
drawn uniformly from a range [—A, A] is 

/A pA 
dxP(x'; t)W(x f -> x) - I dxP(x; t)W(x -> x'). (1) 
A J -A 

Here P{x\t) is the probability of finding the particle at position x at time £, W(x —> x f ) = 
(2A) _1 min (1, exp [— f3U(x f ) + f3U(x)]) is the rate of moving a particle from position x to 
position x\ and /3 = l/(k&T). Eq. ([!]) can be expanded in powers of A (which we assume 
to be small; we also assume that U" • A <C £/'), giving to lowest order 

d t P(x- 1) w -d x (vP(x, t)) + d x (Dd x P(x, t)) . (2) 

This is a Fokker-Planck equation with drift velocity v = — (/3/6)A 2 U f {x) + 0(A 3 ) and dif- 
fusion constant D = A 2 /6 + (9(A 3 ), and corresponds to a Langevin dynamics satisfying an 
Einstein relation —v/D « (3U f {r). We can neglect terms higher order in A provided that 
U(x) changes little in the course of a single move. If this condition holds then Metropolis 
MC moves of single particles in an external potential occur in a dynamically realistic way: 
the drift velocity of the particle is proportional to the force acting upon it, and its diffu- 
sion constant is independent of position. This correspondence also holds in two and three 
dimensions. 

In most situations we are interested in interacting particles, not isolated particles in 
external potentials. Here, too, single-particle Metropolis MC evolution can in many cases 
approximate a realistic dynamics. As an illustration, consider two interacting but otherwise 
isolated particles i and j in one dimension, with positions Xi and Xj. Particles interact 
according to a pairwise potential U(xi — xj). Under the single-particle Metropolis Monte 
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Carlo algorithm described before, the master equation for the separation r = Xi — Xj of these 
particles reads 



/A pA 
dfP(r; t)W(r -> r) - / dfP(r; t) W(r -> r ) 
-A J-A 



(3) 



where f = r' — r and W(r — » r') = (2A) _1 min (1, exp [—/3U(r')+ f3U(r)]). As before, 
expansion of this equation in powers of A yields a Fokker-Planck equation with drift velocity 
v ~ — A 2 /3 £7'(r)/6 and diffusion constant D ~ A 2 /6. These results are proportional to 
those obtained by assuming that i and j are subject to a Brownian motion described by the 
equations 

~/Xi = -d Xt U(r) + r/< (4) 

and 

7^- = -d Xj U(r) +rjj, (5) 

where the noise terms 77^ have zero mean and variance (Vii^Vji^)) = ^k^T^8ij8{t — £'), and 
7 is a friction coefficient. Further, the master equation for the center of mass R = |(o^ + £j) 
of the dimer ij reads 

<%P(i?; t) = [ dfP(R - f/2; i) - f -> z i5 ^) + Xj-r-> Xj)} 

J-A 

- / drP(R; t) [W(xi Xi + f, a^) + Wfa, Xj Xj + r)} , (6) 

J-A 

where W(xi — f x^Xj) = (4A) _1 min (1, exp [— @U(xi — Xj) + f3U(xi — f — Xj)\). Ex- 
pansion of Eq. ^ yields (R) = (since no external forces act on the dimer) and 
(R 2 ) = A 2 /(6 • 4) + 0(A 3 U\r)). Hence in the limit of small displacements A the col- 
lective diffusion constant is independent of the force exterted by i on which is what one 
would conclude by adding Eqs. Q and If, by contrast, trial moves of i and j lead to large 
energy changes, then the dimer diffusion constant (R 2 ) oc J^ A df f 2 exp [— f3U(r f ) + j3U(r)\ 
will be suppressed (potentially strongly so) relative to the overdamped ideal. 

These arguments are difficult to extend analytically to larger numbers of particles, but 
they suggest that single-particle Metropolis MC can mimic a realistic dynamics A) if trial 
moves can be made small enough to render negligible higher-order corrections to particles' 
drift and diffusion coefficients, or B) if collective diffusive modes of motion do not dominate 
a system's behavior. The latter situation might arise if a system is very crowded, and 
particles are able to move only small distances, or if assembly of a structure is dominated by 
the addition and detachment of single particles, rather than by cluster-cluster interactions. 
We illustrate this second scenario in Fig. [jja). 



As suggested by these simple arguments, Metropolis Monte Carlo simulations have been 
shown in several cases to approximate a realistic dynamics. Ref. [H] reports good agree- 
ment between the dynamics of a model protein evolved by integration of Langevin equations, 
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FIG. 1: When explicit collective motion does not (a) and does (b) matter. We show self-assembly 
trajectories for the model protein system introduced in Ref. [T2j . which comprises hard rectangles 
on a two-dimensional substrate (see inset to (a)). Rectangles attract via a rectangular nonspecific 
forcefleld (which encourages fluidlike clustering of monomers) and via directional linkers mediated 
by sticky patches labeled E, S and L (which allow the formation of a square lattice of tetramers). 
Only E-S and L-L pairings are reactive. Here we evolve 1800 monomers at 10% area fraction using 
one of three protocols. These are the virtual- move algorithm described in section |lV| ('collective') , 
with basic displacement scale just larger than the maximum range of particle-particle attractions; 
and the same algorithm in which we reject explicit motion of any cluster larger than a monomer 
('single-particle') or a tetramer ('4-particle'). The lines record the fraction of crystalline particles f c 
(i.e. those with three engaged directional bonds) as a function of the number of simulation timesteps 
£; data are averaged over 8 independent simulations. In (a), the nonspecific attraction is 2.25 k^T 
and E-S and L-L bonds confer 4/cbT of attractive energy. The resulting dynamics comprises 
transient, precritical liquidlike fluctuations that develop crystalline order within them. While large 
structures do move collectively if permitted to do so, crystallization is dominated by detachment 
and attachment of single particles from and to fluidlike blobs, and by rotations of single particles 
within these blobs. As a result, we see little difference between collective- and single-particle 
algorithms: the averages of the respective trajectories lie closer than the trajectories' standard 
deviations. Snapshots, which show about 25% of the simulation box, are taken at 5 x 10 6 timesteps; 
crystalline particles are green. A qualitatively different parameter regime is shown in panel (b). 
Here there is no nonspecific attraction; E-S bonds have a strength of 100 k^T; and L-L bonds have 
a strength of 8.25 k^T. The strong E-S bonds encourage the rapid formation of tetramers and 
fragments thereof, and these clusters subsequently assemble into higher-order structures via L-L 
bonds. We see a large difference between a collective algorithm and a single-particle algorithm, 
because the latter underrepresents tetramer diffusion. Restoring explicit collective moves of only 
4 particles gives rise to behavior similar to that of the fully collective algorithm. Snapshots are 
taken at 7.5 x 10 6 timesteps. 
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and by Metropolis MC using small trial moves. Ref. [13] demonstrates near-quantitative 
agreement between Brownian dynamics and Metropolis MC algorithms used to calculate 
the self-diffusion constant of a colloidal fluid and its crystallization dynamics. Metropolis 
simulations of the crystallization of attractive colloids, in parameter regimes dominated by 
single-particle attachment and detachment, were found to behave like their Brownian coun- 
terparts [H]. Refs. jT5l [16] even found accord between the long-time equilibrium dynamics 
of silica computed with Metropolis MC and with Newtonian (not overdamped) dynamics, 
perhaps because the distinction between diffusive and ballistic dynamics is relatively unim- 
portant in a crowded system. Refinement of the dynamical realism of single-particle MC 
trajectories can be achieved by rescaling the basic MC timestep by acceptance rate | [T3l [T7] , 
or by biasing the choice of trial moves according to the forces acting on particles [18] . 

III. COLLECTIVE MOVES: STATIC LINKING SCHEMES 

There also exist situations in which collective diffusive modes of motion are important 
to a system's evolution, and it is inconvenient to make the Monte Carlo trial step size 
small enough that a single-particle algorithm can mimic realistic collective diffusion. Such 
situations can arise if particles' interactions change rapidly with distance or angle (true of e.g. 
certain model proteins, which might extend several nanometers but interact over distances 
comparable to a nanometer) ; if particles are present at low concentration (because they must 
cross large distances to interact); and if a system's assembly is naturally hierarchical. We 
show in Fig. [j](b) an example in which all three of these conditions holds. 

When faced with such a system, one can restore to a Monte Carlo scheme a degree of 
realistic collective motion by making explicit moves of collections of particles. This is often 
done by recursively linking particle pairs (starting, say, with particles i and j) with a given 
probability, and proposing a collective move of the resulting cluster fT9H23] . This procedure 
connects an initial microstate, /x, with a proposed new one, v. We shall write p^(/x — »> v) for 
the probability of linking i and j in the course of proposing the move \i —> v. 

As Fig. [2] illustrates, particles can be linked statically or dynamically. In a static linking 
scheme (see e.g. [2U [22], EU [25]) the likelihood of linking i and j depends only on their 
properties (e.g. energy or degree of proximity) in the original microstate /i, i.e. p^(/x — » 
v ) = Pijil 1 )- We sketch one example of this procedure in Fig. [2](a). We shall specialize 
our discussion in this section to the simple case of a pairwise potential e^- composed of a 
hard-core part plus an attractive part; one possible static linking scheme is as follows. We 
choose an initial seed particle z, and attempt to link this particle to one of its neighbors j 
with probability 

Pij (n i/) = (n c - ne) lij(fi) (1 - e fte -) , (7) 

where fa is a fictitious reciprocal temperature. The factor 2^(/i) is unity if i and j interact 
(e.g. lie within possible interaction range) in microstate /x, and zero otherwise; the factor 
(n c — rtc) terminates link formation if the number of particles recruited to the cluster, n^, 
exceeds n c , the smallest integer larger than Here £ is a random variable drawn uniformly 
from the interval [0, 1]. If nc exceeds n c then the move is aborted in situ, preventing clusters 
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FIG. 2: Illustration of static and dynamic linking procedures for cluster moves. In (a), nanoparticle 
i is linked to its neighbors j, k and / according to pairwise energies of interaction in the initial 
microstate, p. All recruited neighbors (e.g. j) propose links with their neighbors (e.g. m, n), and so 
on until no particles remain to be tested. In the example shown all particles interact strongly, and 
the entire cluster is chosen to move. The move shown results in proposed new microstate v. In (b), 
i is recursively linked to its environment according to gradients of interaction energies calculated 
by making virtual moves of particles (see text). The proposed move leads to new microstate v. 

from moving with a frequency greater than is physical. 



If the link is succesful then we add j to the moving cluster; if not, we do not, and we do 
not attempt to link i to j again. We continue iteratively to propose links between particles 
in the moving cluster and those with which they interact (as long as we have not tested 
those links before, and provided that those particles are not already members of the moving 
cluster). We stop when we run out of particles to test. We then propose a move (e.g. a 
translation) of the cluster. This defines a new microstate v. To preserve the equilibrium 
distribution it is sufficient to impose the requirement of superdetailed balance [26j , 

p(p)W gen (p is\n)W acc (p -+v\K) = p(v)W gen (is p\TZ)W acc (u -+ p\TZ). (8) 

Here p(p) = e~ l3E ^/Z is the equilibrium weight of the state p (E(p) is the energy of state 
p and Z is the partition function), and W gen (p — >■ v\7l) is the rate of generating a move 
from state p to state v ) given a realization 1Z of links and failed links. This rate contains the 
likelihood of selecting the cluster's displacement or rotation, one factor of Pij(p) for each link 
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formed within the moving cluster, one factor of 1 — Pij{fi) for each link attempted but not 
formed within the cluster, and one factor of 1 —pij(fi) for each link not formed between the 
cluster and its environment. All but the latter set of probabilities equal their counterparts 
for the reverse move. Rearranging Eq. ^ reveals that balance is satisfied by the acceptance 
rate 

W acc (fi i/) = V{C) min (l, e (Pf-P)(E(u)-E^ ? ( 9 ) 

provided that no overlaps occur. If they do, the move is rejected. For infinite fictitious 
temperature, /3f = 0, the likelihood of forming links between the seed i and any other particle 
is zero, and the algorithm executes single-particle moves. For finite values of /3f, collective 
moves can be achieved. We have attached a factor of V(C) < 1 in order to modulate the 
diffusivity of clusters according to their size and shape. 

Such schemes provide a convenient way to effect collective motion in self-assembling sys- 
tems f2l 1251 1271 429] . They allow for precise control of collective motion: we know in advance 
of the move the nature of the cluster to be moved, and we can rotate and translate this 
cluster as desired. Their chief shortcoming, however, is that clusters (and single particles) 
do not move solely according to the potential energy gradients acting on them. Because 
links are conditioned upon energies in the initial microstate, particles interacting strongly 
are likely to be moved in concert, even if relative moves of these particles are favorable. 
The analog of Eq. ^ for a static cluster algorithm yields an effective drift velocity for 
the inter-dimer separation that is not simply proportional to the negative of the potential 
gradient, but is instead proportional to — e^-(r)e /3fe ^ r ). This drift velocity is not consistent 
with a physical dynamics. One suggested consequence of such a bias is sketched in Fig[2](a): 
even though rotation of the upper 6 particles might be desirable (see panel (b)), if particles 
interact strongly then a static algorithm would have trouble forming a cluster of those 6 
particles that does not include the whole of the structure shown. Simulations [30] show 
that proposing relative moves of strongly-attracting particles less frequently than moves of 
weakly-attracting particles can lead to dynamical trajectories substantially different than 
are generated by integrating equations of motion. 



IV. COLLECTIVE MOVES: DYNAMIC LINKING SCHEMES 

By contrast, the idea behind a dynamic cluster-linking scheme [23j is to make a trial move 
of a single particle and to deal iteratively with the consequences of that move. This scheme, 
and certain of its off-lattice generalizations [20l [24], decouple the likelihood of proposing 
relative moves of particles from their interaction energies in the initial microstate, circum- 
venting the chief deficiency of static cluster-linking schemes. In Fig. [2](b) we illustrate one 
possible dynamic cluster-linking scheme [31j, called a 'virtual-move' Monte Carlo algorithm, 
for particles bearing general pairwise interactions . We link particles in a recursive man- 
ner similar to that described above, except that now our linking procedure involves trial 
virtual moves of particles. In detail, we pick a particle i. We choose to form a pre-link of 
particle i and some neighbor j with a probability 



Pijifi -> v) = G (n c - n c )X^(/i)max (0, 1 - <J* u *M-P u i'M) 



(10) 
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that depends on a virtual move (e.g. translation or rotation) of i relative to j. Here Uij(fi) is 
the pairwise energy of the bond ij in microstate /i, and U#j(ii) is the bond energy following 
the virtual move of i (i is returned to its original position following its virtual move). The 
factors 1>ij(n) and B (n c — nc) are as before, and as before the move is aborted if nc exceeds 
n c . Linking particles in this fashion ensures that neighbors exert mutual forces proportional 
to the gradient of their pairwise energies, unlike in static linking schemes. Particle motions 
are correlated at the level of a single move if /3[/^j(/x) — /3Uij(fi) is large. 
We then do as follows. 

• If a pre-link does not form, we label the link ij as unformed, do not add j to the moving 
cluster, and do not consider the link ij again. We then consider another neighbor of i. 

• If the pre-link forms, then 

— we convert the pre-link into a full link with probability /reverse —> v) = 
min ^1, p! J [^y) ) 3 an d add j to the moving cluster, j is then assigned a vir- 
tual move so that it moves with i as a rigid body. To compute /reverse we make 
a reverse virtual move of i (starting from its original position), corresponding to 
the forward virtual move with the sense of rotation or translation reversed. For 
pre-linked particles (indeed, for any two particles internal to the chosen cluster) 



the factor pij(v —> ji) is given by Eq. (10) where i f now refers to i following a 
reverse virtual move. 

— we convert the pre-link into a frustrated link with probability 1 — /reverse (^ ~~ ^ v \ 
In this case j is not added to the moving cluster, and the bond ij is not tested 
again. 

We stop when no more particles remain to be tested, and we move the chosen cluster 
according to the prescribed virtual move, defining a proposed new microstate v. 

To preserve the equilibrium distribution we can balance the rates for forward and reverse 
moves involving a given realization 1Z of 1) internal cluster links, and 2) failed internal links 
that are either unformed or frustrated. By construction of the linking procedure these two 
classes of probabilities cancel from Eq. ([8|. The remaining contribution to Eq. ^ comes 
from unformed links external to the moving cluster, and rearrangement of that equation 
reveals that an appropriate acceptance rate for the collective move is 

W acc (v -+ = V(C) min I 1, [[ e -W^)-u^)) j ? (n) 

provided that no frustrated links lie external to the pseudocluster; the acceptance rate is 
zero if they do. The label (ij) n <+ identifies particle pairs that start (/x) in a noninteracting 
configuration and end {y) with positive energy of interaction (overlapping), or which start 
(/i) in an overlapping configuration and end (y) in a noninteracting one. 

There is no guarantee that such a procedure will result in motion that is dynamically 
realistic in all details. Indeed, there are some features of the algorithm that make precise 
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control of cluster motion impossible. For one, we do not know in advance the nature of the 
moving cluster, making it hard to cleanly separate translational from rotational motion (the 
same is not true of a static linking procedure). For another, collective motion requires that 
both forward and reverse virtual moves of one particle 'recruit' another: if the basic scale 
of virtual displacements is much smaller than particle interaction ranges then such motion 
is unlikely. Correspondingly, if the basic scale of virtual displacements is too large, intra- 
cluster relaxation becomes slow. Some tinkering is needed in order to reach a reasonable 
compromise between the these two processes. 

Even given these difficulties, intuition suggests that by moving clusters according to gra- 
dients of potential energy, and by choosing cluster diffusion constants D(C) in a reasonable 
way, we should reproduce some key features of overdamped motion. (It is also worth noting 
that we have the freedom to scale collective diffusion constants D(C) anisotropically, as is 
physically reasonable, which is not a feature that emerges from simple BD algorithms.) Par- 
ticles should move in a locally realistic fashion and retain some of the collective degrees of 
freedom that single-particle moves ignore. A qualitative comparison between a virtual-move 
algorithm (the version of Ref . J21j , which is the predecessor of the algorithm described here) 
and BD simulations of strongly-attractive discs shows this to be the case, even in circum- 
stances where single-particle moves clearly lack dynamical accuracy [24] . Preservation of 
these important features of real dynamics may be sufficient to determine if the real coun- 
terpart of a model system will assemble well or become kinetically trapped. Testing of a 
virtual-move algorithm against BD simulations of viral capsid self-assembly found that each 
generates similar values of capsid yields for given model parameters [3]]. Since yields depend 
upon both thermodynamics and dynamics, such agreement is encouraging. Other work has 
used virtual-move algorithms (the original version [24J or the one described here) to gener- 
ate dynamical trajectories for self- assembling systems [T21 l32H37| or to thermodynamically 
sample them [38] . 

The cluster schemes described above treat pairwise-interacting particles. However, it is 
possible to use them to effect collective motion of particles bearing multibody potentials [39j , 
which are often encountered in model biomolecules [40j . One way to do so is described in [8] . 
Let's say that the true energy of a system of particles in microstate /i is i?(/i), which may 
contain contributions from multibody potentials. We can nonetheless use the virtual-move 
scheme by assuming that all particles interact via fictitious pairwise potentials (perhaps 
derived from potentials of mean force obtained using particles' true interactions). If we write 
E(n) = {E(n) — U{n)) + U(fi) in the exponentials in the equilibrium weights in Eq. ([8]), 
where U{ji) = \ £V Uij(fi) is the system's total fictitious energy in microstate /x, we find 
the acceptance rate for the virtual-move procedure using the fictitious potentials to be 

W acc (fi -+ v\K) = V(C) min [ 1, e -^ AE - AC/ ) JJ e -mM-UiM) J , (12) 

\ {lj)n++o / 

subject to the same caveats as Eq. (11). Here AE = E(v) — E(/ll) and AU = U{y) — U{n). 
The factor e -£(A£-At/) accounts for the difference between real and fictitious potentials. A 
fictitious linking potential can also be used if the real potential E contains long range inter- 
actions (e.g. 1/r Coulomb interactions) that make direct application of a cluster algorithm 
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inconvenient. In this case, the fictitious potential U could be chosen to account only for the 
short range component of particles' interactions. 

V. CONCLUSIONS 

We have described the use of Monte Carlo algorithms to approximate the overdamped 
dynamics of interacting particle systems. While neither single-particle- nor collective-move 
algorithms are dynamically realistic in all details, recent work shows that they can approxi- 
mate a natural dynamics within a range of model systems. Given the advantages of numerical 
stability and ease of implementation offered by Monte Carlo algorithms over Brownian Dy- 
namics schemes, we suggest that Monte Carlo algorithms can in some cases provide a useful 
and convenient alternative to conventional integration of equations of motion. 
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